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We experimentally observe many-body localization of interacting fermions in a one-dimensional quasi-random optical lattice. We identify the many- 
body localization transition through the relaxation dynamics of an initially-prepared charge density wave. For sufficiently weak disorder the time 
evolution appears ergodic and thermalizing, erasing all remnants of the initial order. In contrast, above a critical disorder strength a significant portion 
of the initial ordering persists, thereby serving as an effective order parameter for localization. The stationary density wave order and the critical 
disorder value show a distinctive dependence on the interaction strength, in agreement with numerical simulations. We connect this dependence to the 
ubiquitous logarithmic growth of entanglement entropy characterizing the generic many-body localized phase. 


Introduction The ergodic hypothesis is one of the central principles 
of statistical physics. In ergodic time evolution of a quantum many-body 
system, local degrees of freedom become fully entangled with the rest of 
the system, leading to an effectively classical hydrodynamic evolution 
of the remaining slow observables ID Hence, ergodicity is responsible 
for the demise of observable quantum correlations in the dynamics of 
large many-body systems and forms the basis for the emergence of local 
thermodynamic equilibrium in isolated quantum systems ElIslHl. It is 
therefore of fundamental interest to investigate how ergodicity breaks 
down and search for alternative, genuinely quantum paradigms in the 
dynamics, and to understand the long-time stationary states that ensue 
in the absence of ergodicity. 

One path to breaking ergodicity is provided by the study of inte- 
grable models, where thermalization is prevented due to the constraints 
imposed on the dynamics by an infinite set of conservation rules. Such 
models have been realized and studied in a number of experiments with 
ultracold atomic gases niiiiTi. However, integrable models represent 
very special and fine-tuned situations, making it difficult to extract gen¬ 


eral underlying principles. 

Theoretical studies over the last decade point to many-body lo¬ 
calization (MBL) in a disordered isolated quantum system as a more 
generic alternative to thermalization dynamics. In his original pa¬ 
per on single-particle localization, Anderson already speculated that 
interacting many-body systems subject to sufficiently strong disorder 
would also fail to thermalize IS). Only recently, however, have con¬ 
vincing theoretical arguments been put forward that Anderson local¬ 
ization remains stable under the addition of moderate interactions, 
even in highly excited many-body states |[^[Tol[TT]. Further theoret¬ 
ical studies have established the many-body localized state as a dis¬ 
tinct dynamical phase of matter that exhibits novel universal behavior 
(I2lllllllll5|[l6l[l7l[l8l[I3|2i^ In particular, the relax¬ 
ation of local observables does not follow the conventional paradigm 
of thermalization and is expected to show explicit breaking of ergodic¬ 
ity. In many ways, the MBL transition is fundamentally different from 
all other known transitions II23I l24l . On one side of the transition er¬ 
godicity prevails and quantum effects decay at long times, whereas on 



Figure 1 : Schematics of the many-body system, initial state and phase-diagram. A. Initial state of our system consisting of a charge density wave, where all 
atoms occupy even sites (e) only. For an interacting many-body system, the evolution of this state over time depends on whether the system is ergodic or not. B. 
Schematic phase diagram for the system: in the ergodic, delocalized phase (white) the initial CDW quickly decays, while it persists for long times in the non-ergodic, 
localized phase (yellow). The striped area indicates the dependence of the transition on the doublon fraction, with the black solid line indicating the case of no 
doublons. The black dash-dotted line represents the experimentally observed transition for a finite doublon fraction, extracted from the data in Fig. 4. The grey arrows 
depict the postulated pattern of renormalization group flows controlling the localization transition. For U = 0, as well as in the limit of infinite U with no doublons 
present ED, the transition is controlled by the non-interacting Aubry-Andre critical point, represented by the unstable grey fixed points. Generically, however, it is 
governed by the MBL critical point, shown in red. C. Schematic showing a visual representation of the three terms in the Aubry-Andre Hamiltonian (Eq. |2|). 
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the other side quantum correlations persist indefinitely. Hence the MBL 
transition sets a sharp boundary between a macroscopic world showing 
quantum phenomena and one governed by classical physics. 

While Anderson localization of non-interacting particles has been 
experimentally observed in a variety of systems, including light scat¬ 
tering from semiconductor powders in 3D II25K photonic lattices in ID 
1^ and 2D ||27l and cold atoms in ID and 3D random (281 EH EOl 
and quasi-random ISD disorder, the interacting case has proven more 
elusive. Initial experiments with interacting systems have focused on 
the superfiuid (3211^ or metal (34) to insulator transition in the ground 
state. Evidence for inhibited macroscopic mass transport was reported 
even at elevated temperatures IMI, but is hard to distinguish from ex¬ 
ponentially slow motion expected from conventional activated transport 
or effects stemming from the inhomogeneity of the cloud. Until now 
conclusive experimental evidence for many-body localization at finite 
energy density has thus been lacking. 

In this paper we report the first experimental observation of ergod- 
icity breaking due to many-body localization. Our experiments are 
performed in a one-dimensional system of ultracold fermions in a bi- 
chromatic, quasi-randomly disordered lattice potential. We identify the 
many-body localized phase by monitoring the time evolution of local 
observables following a quench of system parameters. Specifically, 
we prepare a high-energy initial state with strong charge density wave 
(CDW) order (as shown in Fig. lA) and measure the relaxation of this 
charge density wave in the ensuing unitary evolution. Our main observ¬ 
able is the imbalance X between the respective atom numbers on even 
(Ne) and odd (No) sites 


Ne-No 
iVe + iVo ’ 


( 1 ) 


which directly measures the CDW order. While the initial CDW (X > 
0.9) will quickly relax to zero in the thermalizing case, this is not true in 
a localized system, where ergodicity is broken and the system cannot act 
as its own heat bath (Fig. IB) t35l . Intuitively, if the system is strongly 
localized, all particles will stay close to their original positions during 
time evolution, thus only smearing out the CDW a little. A longer local¬ 
ization length ^ corresponds to more extended states and will lead to a 
lower steady state value of the CDW. The long-time stationary value thus 
effectively serves as an order parameter of the MBL phase and allows us 
to map the phase boundary between the ergodic and non-ergodic phases 
in the parameter space of interaction versus disorder strength. In par¬ 
ticular, in the non-interacting system the CDW vanishes asymptotically 
as oc l/(f^ 1^ . In contrast to previous experiments, which studied the 
effect of disorder on the global expansion dynamics II28II311[^II341[^ . 
the CDW order parameter acts as a purely local probe, directly capturing 
the ergodicity breaking. 

Our system can be described by the one-dimensional fermionic 
Aubry-Andre model Ell with interactions on, given by the Hamil¬ 
tonian 


H = -JY^ (cl^Ci+l,a + h.C.) 

i,a 

+ A cos(27r^i + + U'^ 


( 2 ) 


Here, J is the tunneling matrix element between neighboring lattice sites 
and c\ ^ (Ci^cr) denotes the creation (annihilation) operator for a fermion 
in spin state a G {t, i} on site i. The second term describes the quasi¬ 
random disorder, i.e. the shift of the on-site energy due to an additional 
incommensurate lattice, characterized by the ratio of lattice periodicities 


/3, disorder strength A and phase offset 0. Lastly, U represents the on¬ 
site interaction energy and hi^a = c\ ^Ci^a is the local number operator 
(see Fig. 1C). 



Figure 2: Time evolution of an initial charge-density wave. A charge den¬ 
sity wave, consisting of fermionic atoms occupying only even sites, is allowed 
to evolve in a lattice with an additional quasi-random disorder potential. After 
variable times the imbalance X between atoms on odd and even sites is measured. 
Experimental time traces (circles) and DMRG calculations for a single homoge¬ 
neous tube (lines) are shown for various disorder strengths A. Each experimental 
datapoint denotes the average of six different realizations of the disorder potential 
and the error bars show the standard deviation of the mean. The shaded region 
indicates the time window used to characterise the stationary imbalance in the 
rest of the analysis. 



A/J 

Figure 3 : Stationary values of the imbalance X as a function of disorder A 
for non-interacting atoms. The Aubry-Andre transition is at A/J = 2. Circles 
show the experimental data, along with Exact Diagonalization (ED) calculations 
with (red line) and without (grey line) trap effects. Each experimental data point 
is the average of three different evolution times (13.7r, 17.Ir and 20.5r) and 
four different disorder phases cf), for a total of 12 individual measurements per 
point. To avoid any interaction effects, only a single spin component was used. 
The ED calculations are averaged over similar evolutions times to the experiment 
and 12 different phase realizations. Error bars show the standard deviation of the 
mean. 

This quasi-random model is special in that, for almost all irrational 
[3 (3^, all single particle states become localized at the same critical 
disorder strength NjJ — 2 IITtI . For larger disorder strengths the lo¬ 
calization length decreases monotonically. Such a transition was indeed 
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Figure 4: Stationary imbalance for various interaction and disorder strengths. A: Stationary Imbalance X as a function of interactions U and disorder strength 
A. Moderate interactions reduce the degree of localization compared to the non-interacting or strongly interacting cases. The white dotted lines are contours of equal 
X, while the solid white line is the contour of X matching the Aubry-Andre transition (U = 0 and A/J = 2) extended to the interacting case. It indicates the MBL 
transition. The green dot-dashed line shows the fitted minima of X for each A(36l. Each individual data point (vertices of the pseudo-color plot) is the average of the 
same 12 parameters as in Fig. 3. The color of each square represents the average imbalance of the four points on the corners. All data taken with a doublon fraction 
of ~ 34 (2) %. B: Cuts along four different disorder strengths. The effect of interactions on the localization gives rise to a characteristic ’W’-shape. Solid lines are 
the results of DMRG simulations for a single homogeneous tube. Error bars indicate the standard deviation of the mean. 


observed experimentally in a non-interacting bosonic gas EH. In con¬ 
trast, truly random disorder will lead to single-particle localization in 
one dimension already for arbitrarily small disorder strengths. Previ¬ 
ous numerical work indicates many-body localization in quasi-random 
systems to be similar to that obtained for a truly random potential ESI. 

Experiment We experimentally realize the Aubry-Andre model by 
superimposing on the primary, short lattice (As = 532 nm) a second, 
incommensurate disorder lattice with Xd = 738 nm (thus (3 — Xg/Xd ^ 
0.721) and control J, A and 0 via lattice depths and relative phase be¬ 
tween the two lattices 1^ . The interactions (U) between atoms in the 
two different spin states |t) and \D are tuned via a magnetic Feshbach 
resonance 01 In total, this provides independent control of U, J and 
A and enables us to continuously tune the system from an Anderson 
insulator in the non-interacting case to the MBL regime for interacting 
particles. 

An additional long lattice {Xi = 1064nm = 2As) forms a period- 
two superlattice 13^1^ together with the short lattice and is employed 
during the preparation of the initial CDW state, and during detection 
01 - Deep lattices along the orthogonal directions (A^ = 738nm 
and V± = 36{1)Er), create an array of decoupled ID tubes. Here, 
Er = I (2mAii,) denotes the recoil energy, with h being Planck’s 
constant, m the mass of the atoms and Aiat the respective wavelength of 
the lattice lasers. 

We employ a two component degenerate Fermi gas of atoms, 
consisting of an equal mixture of 25-30 x 10^ atoms in each of the two 
lowest hyperfine states |X, ?tif) = ||?~|) = \i) II’“D — lt)» 
at an initial temperature of 0.24(2) Tp, where Tp is the Fermi tem¬ 
perature. The atoms are initially prepared in a finite temperature band 
insulating state 1401 in the long and orthogonal lattices. We then split 
each lattice site by ramping up the short lattice in a tilted configuration 
1^ and subsequently ramp down the long lattice. This creates a charge 
density wave, where there are no atoms on odd lattice sites but zero, 
one or two atoms on each even site II39II411 . This initial CDW is then 


allowed to evolve for a given time in the 3.0(2)Er deep short lattice 
at a specific interaction strength U in the presence of disorder A. In a 
final step, we detect the number of atoms on even and odd lattice sites 
by employing a band-mapping technique which maps them to different 
bands of the superlattice II41II36I . This allows us to directly measure the 
imbalance X, as defined in Eq. 

Results We track the time evolution of the imbalance X for various in¬ 
teractions U and disorder strengths A (see Fig. 2). At short times the im¬ 
balance exhibits some dynamics consisting of a fast decay followed by a 
few damped oscillations. After a few tunneling times r = /i/(27r J) the 
imbalance approaches a stationary value. In a clean system (A/J = 0) 
and for weak disorder, the stationary value of the imbalance approaches 
zero. For stronger disorder, however, this behaviour changes dramat¬ 
ically and the imbalance attains a non-vanishing stationary value that 
persists for all observation times. Since the imbalance must decay to 
zero on approaching thermal equilibrium at these high energies, the non¬ 
vanishing stationary value of X directly indicates non-ergodic dynamics. 
Deep in the localized phase, where unbiased numerical Density-Matrix 
Renormalisation Group (DMRG) calculations are feasible due to the 
slow entanglement growth, we find the stationary value obtained in the 
simulations to be in very good agreement with the experimental result. 
These simulations were performed for a single homogeneous tube with¬ 
out any trapping potentials 13^ . The stronger damping of oscillations 
observed in the experiment can be attributed to a dephasing caused by 
variations in J between different ID tubes 1411 [^ . 

We experimentally observe an additional very slow decay of X 
on a timescale of several hundred tunneling times for all interaction 
strengths, which we attribute to the fact that our system is not per¬ 
fectly closed due to small losses, technical heating and photon scattering 
|4^[3^. Another potential mechanism for delocalization at long times 
is related to the intrinsic SU(2) spin symmetry in our system 14^ . How¬ 
ever, for the relevant observation times our numerical simulations do not 
indicate the presence of such a thermalization process. 
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To characterize the dependence of the localization transition on U 
and A, we focus on the stationary value of X, plotted in Fig. 3 for non¬ 
interacting atoms and in Fig. 4 for interacting atoms. For non-interacting 
atoms (Fig. 3), the measured imbalance is compatible with extended 
states within the finite, trapped system for A/ J < 2. Above the critical 
point of the homogeneous Aubry-Andre model at A/ J = 2 ITtI . how¬ 
ever, the measured imbalance strongly increases as the single-particle 
eigenstates become more and more localized. The observed transition 
agrees well with our theoretical modeling including the harmonic trap 

EH. 



Time (x) 



U/J 


Figure 5: Growth of entanglement entropy and corresponding slope. A: 

DMRG results of the entanglement entropy growth for various interaction 
strengths and A/J = 5. For long times, logarithmic growth characteristic of 
interacting MBL states is visible. The experimentally used evolution times indi¬ 
cated by the yellow shaded region are found to be in the region of logarithmic 
growth. B: The slope of the logarithmic growth, extracted using linear fits up to 
the longest simulated time (50 r) in A, shows a non-monotonic dependence on 
the interaction strength, which tracks the inverse of the steady state CDW value 
(red line). Error bars reflect different initial starting times for the fit. 

The addition of moderate interactions slightly reduces the degree 
of localization compared to the non-interacting case, i.e. they decrease 
the imbalance X and hence increase the critical value of A necessary to 
cross the delocalization-localization transition (Fig. 4A and B). Impor¬ 
tantly, we find that localization persists for all interaction strengths. For 
a given disorder, the imbalance X decreases up to a value of U ^ 2A 
before increasing again. For large \U\, the system even becomes more 
localized than in the non-interacting case. This can be understood qual¬ 
itatively by considering an initial state consisting purely of empty sites 
and sites with two atoms (doublons): for sufficiently strong interactions, 
isolated doublons represent stable quasiparticles as the two atoms cannot 
separate and hence only tunnel with an effective second-order tunneling 


rate of Jd = ^ J 1^1^ . This strongly increases the effective 

disorder oc Jd ^ A/J and promotes localization. In the experi¬ 
ment, the initial doublon fraction is well below one EH and the density 
is finite, such that we observe a weaker effect. We find the localization 
dynamics and the resulting stationary values to be symmetric around 
f/ = 0, highlighting the dynamical U ^ —U symmetry of the Hubbard 
Hamiltonian for initially localized atoms 1461 . The effect of interactions 
can be seen in the contour lines (Fig. 4A, dotted white lines) as well 
as directly in the characteristic ‘W’ shape of the imbalance at constant 
disorder (Fig. 4B), demonstrating the re-entrant behaviour of the sys¬ 
tem l22l . This behaviour extends to our best estimate of the localization 
transition, which is shown in Fig. 4A as the solid white line. 



U/J 

Figure 6: Stationary imbalance X as a function of interaction strength dur¬ 
ing loading. Data taken with disorder A/J = 3. The loading interactions of 
<^ioad = —89(2)ao (attractive, where ao denotes Bohr’s radius), 0(1) ao (non¬ 
interacting) and 142(1) ao (repulsive) correspond to initial doublon fractions of 
51(1)%, 43(2)%, and 8(6)%, respectively (21- Each X value is the average of the 
same 12 parameters as in Fig. 3. Error bars show the standard deviation of the 
mean. Solid lines are guides to the eye. The grey shaded area spans the limiting 
cases of 0 and 50% doublons, simulated using DMRG for a single homogeneous 
tube. 

We can gain additional insight into how localization changes with 
interaction strength by computing the growth of the entanglement en¬ 
tropy between the two halves of the system during the dynamics, as 
shown in Fig. 5A. For long times, we observe a logarithmic growth of 
the entanglement entropy with time as S{t) = A^offset + s* ln(t/T), 
which is characteristic of the MBL phase ifT^ [T^ . The slope s* is 
proportional to the bare localization length which in a weakly in¬ 
teracting system in the localized phase corresponds to the single particle 
localization length. In general, is the characteristic length over which 
the effective interactions between the conserved local densities decay 
(TTlITSl and connects to the many-body localization length ^ deep in the 
localized phase. In contrast to however, * is expected to remain finite 
at the transition 1^ . We find s* to exhibit a broad maximum for inter¬ 
mediate interaction strengths (Fig. 5B), corresponding to a maximum 
in the thus inferred localization length. This maximum in turn leads to 
a minimum in the CDW value. The characteristic ‘W’ shape in the im¬ 
balance is thus directly connected to the maximum in the entanglement 
entropy slope, as both are consequences of the maximum in localization 
length. Equivalent information on the localization properties as obtained 
from the entanglement entropy can be gained in experiments by moni¬ 
toring the temporal decay of fiuctuations around the stationary value of 
the CDW (3^. While we do not have sufficient sensitivity to measure 
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these fluctuations in the current experiment, we expect them to be ac¬ 
cessible to experiments with single site resolution 147)1481 . 

To systematically study the effect of the initial energy density on the 
MBL phase, we load the lattice using either attractive, vanishing or re¬ 
pulsive interactions (Fig. 6), thereby predominantly changing the num¬ 
ber of doublons in the initial state 1^ . Since the initial state consists of 
fully localized particles only, the local energy density is directly given 
by the product of interaction strength U and doublon density. We And 
that for an interaction strength during the evolution of \U/ J\ < 6 the 
energy density does not signiflcantly affect the localization properties, 
proving that MBL persists over a wide energy range. For \U/ J\ > 8, 
localization properties depend signiflcantly on the doublon fraction be¬ 
cause of the second emerging energy scale Jd, sls discussed above. For 
the case of repulsive loading, i.e. a low fraction of doubly occupied sites, 
the imbalance for f// J = 0 and strong interactions match within error. 
Indeed, a rigorous mapping can be made between the non-interacting 
system and the dynamics in the doublon free subspace at strong inter¬ 
actions \U/J\ oo (361 . At very large interactions and high dou¬ 
blon fractions, the additional long timescales start to also compete with 
heating and loss processes, rendering the deflnition of stationary states 
challenging. 

Conclusion We have created and characterized the many-body local¬ 
ized phase of a system of interacting ultracold fermions in a quasi¬ 
random optical lattice. The dependence of the MBL phase-transition 
point on interactions was measured by studying the evolution of an arti- 
flcially created charge density wave. Moderate interactions have a delo¬ 
calizing effect and increase the critical disorder strength. We have also 
shown the MBL phase to be stable over a wide range of energy densities. 

Our experimental demonstration of ergodicity breaking due to 
many-body localization paves the way for many further investigations. 
An interesting extension would be to use ‘true’ random disorder created 
by e.g. an optical speckle pattern, as has been used to study Anderson 
localization Another important next step is extending the present 
study to higher dimensions. Additional insight can also be gained by 
analyzing the full relaxation dynamics of local observables 
in an experimental setup featuring single site resolution (4711^ . For 
instance, the decay of fluctuations of X with time could be directly mea¬ 
sured, providing an even more direct connection to the entanglement 
entropy. Another important direction for future investigation is the ef¬ 
fect of opening the system in a controlled way. This could be done, 
e.g. by adding a near-resonant laser to deliberately enhance photon scat¬ 
tering or by employing a Bose-Fermi mixture, in which excitations of 
the BEC form a well controlled bath for the fermions. This will allow 
a systematic study of the critical dynamics associated with the MBL 
phase transition, where the bath relaxation time now provides the only 
scale. Such a study would also allow the MBL phase to be clearly dis¬ 
tinguished from classical glassy dynamics. The latter, unlike MBL, is 
insensitive to coupling of the system to an external bath. 
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Supplementary Material: 

Supplementary Text 

General sequence The experimental sequence for producing degener¬ 
ate Fermi gases of via sympathetic cooling with in a magnetic 
quadrupole and optical dipole trap followed by evaporative cooling has 
been described elsewhere ITOl . Final clouds contain typically 50-60k 

atoms at a temperature of 0.24(2) times the Fermi temperature Tf, 
with 50(3)% of atoms in each of the |t) and |T) spin states. We tune 
the interactions between the two spin states via a Feshbach resonance at 
202.IG |49l- All experiments described in this work employ absorption 
imaging after 8 ms time of flight (TOF) along the y axis {z axis is ori¬ 
ented along the vertical direction). This is orthogonal to the x axis along 
which the superlattice is aligned and all experiments take place. For the 
experiments shown in Fig. 3, where we employ a single spin state, we 
set the magnetic held early in the evaporation sequence to 469 G, which 
is the centre of a Feshbach resonance between |t) atoms and 
atoms in the \F — 2^mF — 2) state ll50l . This results in all |t) atoms 
being lost (along with some ®^Rb atoms). By adjusting the cooling se¬ 
quence slightly we are able to prepare a sample with 98(5)% of atoms in 
the IT) state with the same temperature and atom number as in the usual 
sequence. 

Lattice Details: 

Lattice setup All lattice potentials result from standing waves cre¬ 
ated by retro-reflected single-frequency laser beams along ihQx.y and z 
axes. All beams have a Gaussian profile and are focussed down to 1/e^ 
waists of ~ 150 ym at the positions of the atoms. As our cloud has a 
FWHM of ~ 50 ym in the horizontal plane and ~ 15 ym in the vertical 
direction, this results in slight variations of the lattice parameters J, U 
and A across the cloud (see below). 

Superlattice lock In our superlattice setup the 1064 nm laser 
(long lattice) serves as a master oscillator and its frequency is locked 
to a thermally stabilized and acoustically isolated Fabry-Perot cavity. 
The absolute short term stability of the master oscillator is 65 kHz over 
100 ms and is characterized through the residual locking error. The short 
lattice laser at 532 nm is offset locked relative to the second harmonic 
of the long lattice laser. The offset lock frequency sets the relative phase 
of the superlattice potential 1^ . We obtain a relative line width of 
the offset lock of 1.6 MHz, which corresponds to a phase fluctuation of 
(50 = 14mrad. This also dominates the absolute short term stability of 
the short lattice laser. 

Phase adjustment of the superlattice In order to calibrate the 
phase of the superlattice, we prepare a spinpolarized band insulating 
cloud in a combination of long and perpendicular lattices and subse¬ 
quently split each well non-adiabatically by adding the short lattice. 
Time of flight pictures of such arrays of double wells show the typical 
double-slit interference pattern only for phase settings of the superlat¬ 
tice that result in a symmetric splitting. By determining the difference 
in offset frequencies necessary to reach neighboring symmetric conflg- 
urations, we can precisely calibrate the phase 0 of the superlattice. 

Disorder lattice We use a Coherent MBR-110 Ti:Sa laser oper¬ 
ated at 738 nm to create the disorder potential and the perpendicular 
lattices. This laser is stabilized via its internal reference cavity, which 
can be externally tuned to change the laser frequency. An additional ex¬ 
ternal temperature stabilized and vibration isolated Fabry-Perot cavity 


is used to measure the linewidth and to calibrate the frequency changes 
necessary to set the phase 0 of the disorder lattice. The phase 0 is al¬ 
ways set before the lattice sequence starts and kept constant throughout 
the experiment. 

Note that, strictly speaking, the incommensurate ratio 0 532/738 

should be an irrational Diophantine number EB for there to be a sharp 
localization transition in the Aubry-Andre model. However, for a real¬ 
istic flnite system, it is sufficient for the actual period of the combined 
disorder and regular lattices to exceed the system size ED. 

Disorder strength The disorder strength A in the tight-binding 
limit depends on the lattice depth of the primary lattice Vs — Ss •F^i?, 532 , 
the disorder lattice depth Vd = Sd ■ Erjss and the wavelength ratio 
0 532/738. It is given by |52l 

A — ksER^532 
~0.2 • Sd • Er^532 


I 


dx cos{2/3ksX \ w {ksx)f‘ 


(3) 


where w (x) denotes the Wannier function of the unperturbed short lat¬ 
tice, ks = 27r/As, and Ss and Sd are the depths in terms of Er of 
the short and disorder lattices, respectively. Here, A is given in units 
of Er ^332 = h X 17.64 kHz. We use a primary lattice depth of Ss = 8 
yielding an undisturbed tunnel coupling of J = /i x 540 Hz. Combining 
these numbers, we obtain a conversion factor for 


A 
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0.2 • Sd • Er^332 
/i-550 Hz 


6.67 • Sd- 


(4) 


Long-Term evolution For extremely long hold times (t ^ 100 ms 
or 500t) we observe an additional decay of the imbalance, along with 
the total particle number. This decay can be attributed to unwanted ad¬ 
ditional effects, such as off-resonant light scattering and collisions with 
background gas atoms, as well as technical imperfections including laser 
noise, frequency instabilities or vibrations and other pointing instabili¬ 
ties. All of these effects violate the condition of a perfectly closed sys¬ 
tem and will therefore lead to a vanishing imbalance in the long-time 
limit. Since the relevant timescale is, however, more than an order of 
magnitude larger than the longest times used in the experiment, these 
effects have a negligible impact on the reported observations. 


Trap effects on non-interacting dynamics Compared to the homo¬ 
geneous Aubry-Andre model, the experimental results are modifled due 
to the presence of the 3D harmonic trap. All experimental results are au¬ 
tomatically averaged over many inequivalent ID tubes and, furthermore, 
the atoms in each tube experience a conflning harmonic potential. 

In the presence of a longitudinal harmonic trapping potential (trap¬ 
ping frequency 27r x 60 Hz), the sharp transition at A/J = 2 

expected in the homogeneous case oil becomes smeared out (see Fig. 
3), as even in a weak trap the localization length is bounded from above 
by the flnite cloud size. In a sufficiently strong trap, particles situated 
far away from the trap center will experience offset energies that are 
bigger than 4 J and thereby become localized to their respective side of 
the trap, even in the absence of any disorder. Nonetheless, although the 
trap causes some degree of localization even in the absence of disor¬ 
der, we still observe a clear crossover from trap dominated behavior to 
disorder-induced localization. 

The flnite-size of the used lattice beams (waists 150 ym) re¬ 
sults in an additional small but flnite variation of J between the central 
and off-center tubes which increases with the distance from the beam 
center. More than 90% of the atoms see a less than 20% change in the 


7 




Lattice Loading )>j^CDW Preparatior^^Evol.^)> Imbalance Detection 



Figure 7: Full experimental sequence. Schematic showing lattice depths, phase and Feshbach field ramps for loading, CDW preparation, evolution and detection 
of the imbalance. See text for details. 


hopping term. The Rayleigh length of the lattice beams is long enough 
that there is no appreciable change in tunneling along the tubes. The 
above variation in J mainly gives rise to a dephasing between the dy¬ 
namics in individual tubes, and is responsible for the more pronounced 
damping of imbalance oscillations compared to the homogeneous case 
(see Fig.[^. In the ideal homogeneous system with zero disorder and 
no variation in J, the imbalance would decay according to a Bessel func¬ 
tion 15^ . 

Lastly, since the quasi-random potential is produced by superimpos¬ 
ing a disorder lattice on top of the primary lattice, this additional poten¬ 
tial gives rise to a longitudinal modulation of the hopping rate. This 
can be modeled as a sinusoidal variation of the same periodicity as the 
quasi-random potential (i.e. as sin(27r/3z + 0 + 0diff), where 0diff is the 
phase difference between the modulation of on-site energy and tunnel¬ 
ing. For A/J = 3, the amplitude of this modulation is about 8% of the 
hopping term. This modulation does not change the average hopping in 
a tube and hence does not change the average imbalance. However, it 
further increases the damping of the imbalance oscillations. 


CDW preparation sequence The preparation sequence is illustrated 
in Fig.|7] To prepare the initial CDW state we start with a 3D band in¬ 
sulating state in a 20 Er deep long lattice and 36 Er deep orthogonal 
lattices. The short lattice is then ramped up in 200 jis to 20 Er with 
the superlattice phase set to 0 ~ 37r/4. This creates an array of tilted 
double wells, where the tilt offsets are large enough for all atoms to 
be loaded into the even wells. This tilt additionally suppresses all tun¬ 
nelling while we set the interaction strength U by ramping the Feshbach 
field to the desired value in 25 ms. The long lattice is then switched off 
in 10 /is while the disorder lattice is simultaneously ramped up. Finally, 
the short lattice is ramped down to SEr in 10 /xs, thereby enabling tun¬ 
neling and initiating the dynamics. The transverse lattices remain deep 
to decouple individual tubes. 

This preparation sequence yields 96(2)% of atoms on even sites in 
the short lattice. We attribute the small residual number of atoms sit¬ 
ting on odd sites to non-adiabaticities in the splitting process, tunneling 
events during the preparation sequence, higher band occupancies and 
detection imperfections. 






























































Side-resolved detection sequence The detection sequence to deter¬ 
mine the number of atoms on even and odd lattice sites (and hence T) 
is shown in Fig. [T] Once the evolution in the 8 Er short lattice with 
disorder is complete, we increase the short lattice in 10 /xs to 20 Er to 
minimize the tunneling rate. The long lattice is then ramped up to 20 Er 
in 10 /is, while the disorder lattice is simultaneously ramped down. With 
the spatial distribution of atoms now frozen in the tilted superlattice, we 
ramp the Feshbach field to the non-interacting (U = 0) point in 21 ms, 
to ensure that there are no unwanted interaction effects during the subse¬ 
quent band-mapping and imaging process. Next, the long lattice depth 
is further increased to 88 Er in 20 /xs, transferring atoms from the sec¬ 
ond band into the third band of the superlattice via an avoided crossing 
E). Any atoms on odd sites of the short lattice (the higher side of the 
tilted superlattice double wells) are now in the 3rd band, while atoms 
from even sites of the short lattice remain in the lowest band of the su¬ 
perlattice. After 10 /xs of hold time, the short lattice is ramped to 0 Er 
in 200 /xs, which is slow enough to be adiabatic with respect to inter¬ 
band transitions. Since the ordering of the bands remains unchanged 
during this ramp, atoms originally on even (odd) sites in the short lattice 
end up in the 1st (3rd) band of the long lattice. Finally, we implement a 
bandmapping sequence followed by 8 ms TOF and absorption imaging, 
from which we can count the number of atoms in each band directly, 
providing us with Ne and No. A few atoms remaining in the second 
band due to an imperfect transfer are counted together with the third 
band. 


be converted into Feshbach molecules El. Since the molecular tran¬ 
sition is out of resonance of the imaging light, we can determine the 
doublon fraction directly by comparing the number of atoms detected 
at the two different Feshbach ramp endpoints. The number of atoms 
is measured using standard absorbtion imaging after 6ms TOF. Fig. 
shows the doublon fraction plotted against the loading interaction aioad 
during the lattice ramp up. Error bars show the standard deviation. We 
quote the loading interactions in units of ao rather than U/J, as the lat¬ 
tice configuration at the end of the loading differs from the one during 
the evolution. In the lattice configuration used for time evolution of the 
CDW, a scattering length of a = 7.105 ao corresponds to f// J ~ 1 . 

We also investigate the effect of the initial temperature on the dou¬ 
blon fraction. The initial temperature is changed by reducing the dura¬ 
tion of the final sympathetic cooling stage and increasing the final dipole 
evaporation depth. This has the effect of making the cooling less effi¬ 
cient while keeping the final atom number relatively constant. Fig. 

shows the doublon fraction plotted as a function of the temperature of 
the initial state, expressed in units of the Fermi temperature Tp, for three 
different loading interactions. In the case of loading with attractive in¬ 
teractions or in the non-interacting case, we observe higher temperatures 
to reduce the resulting doublon fraction, since the increased entropy and 
average kinetic energy per particle reduce the average density and render 
the formation of doublons less favourable. For atoms which are loaded 
repulsively, there are very few doublons for any temperature and thus no 
real dependence is observable within our error bars. 
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Figure 8: Measured fraction of atoms on doubly occupied sites as a function 
of interaction strength during loading. The initial state prior to loading into 
the lattice has a temperature of T/Tp = 0.24(2). Error bars show the standard 
deviation. 


Figure 9: Measured fraction of atoms on doubly occupied sites as a function 
of initial reduced temperature for three different loading interactions. Error 
bars show the standard deviation. 


Doublon fractions In order to characterize our initial state, we mea¬ 
sure the fraction of atoms forming doubly occupied sites, termed the 
doublon fraction, by converting doublons into molecules by crossing 
the Feshbach resonance fdO). Experimentally this involves loading into 
a 20Er deep long lattice and 36 Er deep orthogonal lattices as per the 
standard sequence with loading interactions given by aioad, then ramping 
the long lattice to 40(1)Er in 200 /xs to supress hopping in all three di¬ 
mensions, which preserves the atomic distribution for the duration of the 
doublon detection sequence. The Eeshbach field is then ramped down 
to either 206.3G (= -117(2)ao, still above the resonance) or 198.3G on 
the molecular side of the resonance. When the Eeshbach field ramp 
crosses the resonance, any pairs of atoms on the same lattice site will 


Imbalance vs energy density We investigate the effects of both load¬ 
ing interactions and initial temperature on the imbalance X in the many- 
body localized regime at A/ J = 2.3(4), close to the localization tran¬ 
sition. While changing the interactions during loading predominantly 
affects the number of doublons for a given temperature fSS], changing 
the temperature of the initial state affects both doublon fraction and av¬ 
erage density. Eig. shows the effect of loading interactions on the 
imbalance X for three different interactions during evolution. 

Eor vanishing or weak interactions during evolution, we observe no 
infiuence of loading interactions, highlighting the robustness of many- 
body localization in a lattice with respect to energy density. However, 
for strong repulsive interactions during evolution, we see a significant 
effect of the loading interactions on the imbalance X. Increasing the 
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number of doublons by loading attractively increases the imbalance, 
since the additional energy scale Aj Jd of the doublons brings the sys¬ 
tem much deeper into the localized phase. In contrast, for repulsive 
loading, doublons are strongly suppressed. In the limit of strong in¬ 
teractions (during evolution) and the absence of doublons, interacting 
fermions can be exactly mapped to spinless (non-interacting) fermions 
(see below). This is precisely what we observe in our experiment for 
strong repulsive loading interactions and thus vanishing doublon frac¬ 
tion: the imbalance observed for strong interactions during evolution 
converges towards the case of a non-interacting evolution. 



Scattering Length During Loading (a^) 


Figure 10: Measured steady state imbalance as a function of loading inter¬ 
actions. Data taken with A/J = 2.3(4) at a temperature of TjTp = 0.24(2) 
prior to loading. Each X value is the average of the same 12 parameters as in Fig. 
3. Error bars show the standard deviation of the mean. 

For intermediate interactions Vj J 5, we cannot detect any no¬ 
table variation of the imbalance X within errors for all loading inter¬ 
actions. Additionally, in Fig. we investigate the effect of loading 
temperatures on X. Fig. Elh shows the case of attractive loading at 
dioad = -89(2) ao, while Fig. shows the case of repulsive load¬ 
ing interactions aioad = 101(1) ao. In all cases, increased temperatures 
reduce the effect of interactions due to the decreasing average density, 
connecting the interacting cases asymptotically to the behaviour of non¬ 
interacting fermions. In the case of attractive loading, higher tempera¬ 
tures in addition reduce the number of doublons, which for strong repul¬ 
sive interactions during the evolution reduces the localization. For inter¬ 
mediate interactions, the effect of initial temperature is less pronounced, 
but for both loading interactions increasing initial temperatures connect 
the minimum of localization (Fig. 6) to the case of no interactions , see 
Fig.[n] 


Fitting of the imbalance vs interaction curves We extract the min¬ 
ima of the imbalance as a function of interactions U/J by employing 
the following phenomenological fitting function: 


X{U) = \ci\e 



|c2|e 


""2 / -p offset. 


(5) 


Here ci, C 2 , wi, W 2 and offset are free fitting parameters. We have 
implemented the U ^ —U symmetry to the function by fixing the cen¬ 
ter of both Gaussians to zero. The positive amplitude Gaussian with 
a smaller width captures the reduction of the imbalance up to the mini¬ 
mum, while the decreasing contribution of the negative amplitude Gaus¬ 
sian resembles the increase for strong interactions. The minima of this 
function are used to generate the green dot-dashed lines in Fig. 4A. 


Time evolution of the CDW in a non-interacting model We con¬ 
sider an initial state with density wave order, i.e. 

^0 = z = {Ne - No)J{Ne + No) ^ 0, (6) 

I 

where L is the total number of sites, Ne (No) is the number of fermions 
on even (odd) sites, and . .)q represents a trace over the initial density 
matrix. This state is evolved with a (spinless) non-interacting fermionic 
Hamiltonian 


H = - J + h.c.) + 

i i 


(7) 


The disorder Vi could represent true randomness or a quasi-periodic po¬ 
tential as in the experiment. To obtain the density wave order at later 
times we need to compute the expectation values of the local densities 
at later times, i.e. ^ — (hi(t))^. We can write the density 

operators fii (t) explicitly as 

ni{t) = cl{t)ci{t) = 'Y^o{3,kt)Go{j' ,l\t)c\cy, (8) 
33' 


where Go(j, l]t) — (J I I 0 the non-interacting propagator 

from point j to point I and XL is the first-quantized Hamiltonian. Tak¬ 
ing the trace over the density matrix leaves only the terms with j — j' 
yielding 

(n;(^))o = ^0 |Go(2m,, (9) 

m 

where we have assumed an initial state with particles only on even sites, 
i.e. (n 2 m+i)o = 0- 

Therefore the CDW or (imbalance) at time t is 

X(i) = 3y^(-l)*|Go(2m,Z;^)|^ (10) 

l,m 

which can be calculated numerically 


Details of system modelling via exact diagonalisation In order to 
calculate the red line in Fig. 3, we employ Eq. to calculate the time 
evolution of the initial CDW in a single tube in the presence of the har¬ 
monic trap. The initial CDW is modeled as a collection of randomly dis¬ 
tributed localized particles, whose distribution follows a Gaussian with 
(gx — 21 /xm) and matches the initial imbalance of Xq — 0.96(2) ob¬ 
served in the experiment. We extract the stationary imbalance by aver¬ 
aging the imbalance in the evolution time interval from 15 r to 30 r. In 
addition, the calculations are averaged over twelve disorder realizations 
corresponding to twelve equally spaced phases 0 G [0, 27r] in Eq. §. 

In order to correctly model the dynamic evolution shown in Eig.p^ 
we additionally integrate the results over many individual ID tubes with 
tube-dependent hoppings in order to fully model the experimental situ¬ 
ation. The widths and amplitudes of the corresponding Gaussian distri¬ 
butions were chosen in accordance with the experimental cloud sizes of 
cFx,y = 21 fjm and Gz — 5.5 /xm that were obtained via in-situ imaging. 
We included both the hopping perturbation resulting from the disorder 
potential Vi — A cos(27r/3z + 0) as well as the effect of the finite waist 
of the lattice beams. The strength and phase difference 4>di //of the dis¬ 
order induced hopping modulation are extracted from fits to the exact 
position dependent hopping rates. 
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Figure 11: Measured imbalance as a function of initial temperature. Data shown for A attractive loading and B repulsive loading. Each X value is the average 
of the same 12 parameters as in Fig. 3. X error bars show the standard deviation, while Y error bars show the standard deviation of the mean. 


Stationary imbalance value in the localized phase We can give an 
analytic estimate for the steady state value at long times in the local¬ 
ized regime. Time averaging the wave function of a particle originat¬ 
ing from a single site should lead to a form which is exponentially 
decaying over a localization length Hence, we can approximate 
~ I + where the prefactor 1+ 6"’"*' 

ensures that the initial site of the wave function is an even site. Note that 
this form holds for ^ ^ d, i.e. a localization length much larger than a 
lattice constant d = As/2. Inserting this expression into Eq. (To| for the 
density wave order yields 


I(ooj ~ - 


IV 


— iTT(l — V) I 
e + e 






— inr I —i^r —ini 
e + e 2 g 




r 




( 11 ) 


where in the last equality we changed to a sum over relative r and center 
of mass coordinates 1. The second term vanishes and we find (for large 
^'>d) 


X(oo) 


1 

L 



\r\/^ 




( 12 ) 



Figure 12: Time evolution of a non-interacting initial charge-density wave. 

We experimentally measure (circles) the time evolution for various disorder 
strengths and compare it to an exact diagonalisation simulation, taking into ac¬ 
count the trapping potential and the variations in tunneling. Experimental error 
bars show the standard deviation of the mean, while the shaded region indicates 
the standard deviation of simulation runs for various disorder phases. 


This formula establishes the connection between the CDW station¬ 
ary value and the localization length close to the transition in a non¬ 
interacting system. 

Mapping from hard-core to non-interacting fermions We consider 
the dynamics of the (spin-1/2) Hubbard model in the limit of infinite 
(on-site) U when the initial state has no doubly-occupied sites. In this 
case, the time evolution can be mapped to that of non-interacting spin¬ 
less fermions. 

In the limit U/J (X) the Hilbert sub-space with no doubly- 
occupied sites is closed under unitary time evolution. Double occu¬ 
pancies cannot be generated in the dynamics while conserving energy 
because they incur an infinite energy cost. Hence, we can replace the in¬ 
teraction with a hard-core constraint implemented by the operator Pq, 
which projects out the doubly-occupied subspace, 

H = [PG^iA,i+XG + h.c.] + Vihi. (13) 

i(j i 

Let us denote the ordering of the fermion spins from left to right on the 
chain with {cr} = {cri,cr 2 ,...,crAf}, where N is the total fermion num¬ 
ber. {cr} is an invariant of the dynamics because in order to exchange the 
spins of two particles with opposite spin orientation, the particles need 
to meet on the same site, which is prohibited by the hard-core constraint. 
Now, suppose we start the time evolution with an initial state having a 
definite spin ordering {a}. Since this is an invariant of the dynamics, 
we can describe the time evolution within this sector with a Hamiltonian 
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H[{a}] in which the spin ordering enters as a parameter. The Hilbert 
space in which this Hamiltonian acts can be faithfully represented by 
the Fock space of spinless fermions | ni,..., because given 

the fermion positions their spins are uniquely determined by the con¬ 
figuration {cr} which parameterizes the space. Hence we can write the 
Hamiltonian which acts within the {a} subspace as a Hamiltonian of 
spinless fermions 

H[W}] = h-c-] + E fi ■ (14) 

ij i 

Note that the projection of doubly-occupied sites is automatically taken 
into account using the representation with spinless fermions. 

The mapping of the time evolution within a given spin-ordering sec¬ 
tor {cr} can also be understood as a result of a non-local unitary change 
of basis 

L 

U[{(y}] = H exp , (15) 

j=i 


where N^j denotes the particle number to the left of point j, — 


V’{<T} [rij] = 


if an 
if an 


(16) 


In other words, the transformation counts the number of particles to the 
left of point j to find its order in the chain rij . It then rotates the spin by 
TT around the x axis if in the configuration {cr} the nj-th spin is The 
transformed fermions are then all t and we can drop the spin index. 

The dynamics is thus described by the same Hamiltonian of non¬ 
interacting spinless fermions regardless of the initial ordering of the 
spins {cr}. In particular, the time dependence of the CDW order can 
be generated by unitary evolution of a spinless non-interacting fermion 
model and is independent of the initial spin ordering. 


Time evolution of the CDW in an interacting model We consider 
the following ensemble of many-body states as initial states for the nu¬ 
merical calculation: 

l^0>=^(4n^t) (cLj (17) 

m 

where \Q) describes the vacuum and riia = 0,1. In order to gener¬ 
ate the ensemble of initial states, the occupations riia = 0,1 are chosen 
from an appropriate joint probability distribution p{n ^, n^) such that the 
average site occupation and the average doublon occupation are fixed. 
This state is time evolved with the interacting Hamiltonian 

H =-J E(4.aCi+l,a + h-C.) + E E (1^) 

i,cr i 

using the time-evolving block decimation (TEBD) 15^ [57l on matrix 
product states implemented with the ITensor Library 1581 . 

For our numerical simulations we use a quasi-periodic potential 
Vi = Acos(27r/3i + 0) with /3 0.721 as in the experiment. How¬ 

ever, we do not take into account the periodic trap, but rather consider 
a system with L = 40 sites and open boundary conditions. When sim¬ 
ulating the time evolution for a given set of system parameters, i.e. A, 
U, site occupation, and doublon density, we average over both the ini¬ 
tial state and the disorder realization (through sampling of the random 
phase 0). The total number of runs for a given set of parameters is 
100 = 10 X 10, i.e. about 10 initial states for every disorder real¬ 
ization. For the singular value decomposition, we use a cutoff of 10“^ 
and allow for a maximum bond dimension of mmax = 100. Using a 
time step of At = 0.1 [h/J], this leads to a maximum truncation error 
^ 10“^- 


Fluctuations of the imbalance In the main text we presented TEBD 
(i.e. DMRG) calculations showing logarithmic growth of the entangle¬ 
ment entropy with time, which is characteristic of the many-body local¬ 
ized state. The entanglement entropy, however, is a non-local quantity 
and therefore not directly measurable. In this section we show that the 
same information about the many-body localized state can be inferred 
from the temporal fluctuations of the measured imbalance X (CDW or¬ 
der parameter) around its stationary value. 

We focus on the system deep in the localized phase, where it can 
be described using an effective model of local integrals of motion. Eor 
simplicity we initialize the system in a CDW state with exactly one par¬ 
ticle on every even site. We now partition the lattice into two-site dimers 
and, as a starting point, take the dimers to be disconnected. We define 
O'! = — fi 2 i+i and af = which within the 

accessible Hilbert space act like Pauli spin operators. In terms of these 
operators the imbalance (or CDW order) reads 

(19) 

i 

Eor a system with strongly localized particles (i.e. (f* ^ 1) we will 
continue to use an effective model of decoupled dimers, where from 
now on L denotes the number of dimers. We take into account the tail of 
the single particle wave-function, which extends outside of the dimers, 
only through the interaction it induces between the local spin operators. 
Hence an effective model is 

H = '^hi • ai Vija-d-j, (20) 

where hi = Jx -\- CiZ, af = {hi ■ ai)/\hi\ are the integrals of motion 
and the interaction is given by Vij ~ U exp (—— Xj\/^^). 

Eor the non-interacting case, only the first term of the effective 
model is non-zero, and we directly find the expectation value of the spin 
operator on site i 


{af (t)) = cos^ Oi + sin^ Oi cos{uit), (21) 

where uji = and cos^ 0 = ef/ujf ^ 1 — (J/a)'^. Thus, 

we find the time-averaged stationary value of the imbalance is X = 
Z temporal fiuctuations around this average are 

{SX{t)^)j. = / ^ ^ sin^ Oi sin^ Oj cos{uJit) cos{ujjt) \ 

^ ^ ( 22 ) 

where (.. .)^ represents averaging over a long time window T. We con¬ 
clude that without interactions the fluctuations of the imbalance remain 
on average constant in time with (5Xrms = ((5X(t)2)^ ~ 1 /\/L. 

When interactions are present the local contribution to the density 
wave (erf (t)) is affected by the time evolution of the other spins. The 
pseudo spin af gradually becomes entangled with a growing number 
of other spins. Specifically, because the interactions Vij decay exponen¬ 
tially with distance, the number of spins that become significantly entan¬ 
gled with a given spin at time t is l{t) ~ S{t) = SoEset + s* \xi{Ut/h). 
We will see below how this affects the fluctuations of {af) and ulti¬ 
mately contributes to the temporal fiuctuations of the order parameter 
X. 
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Let us calculate the reduced density matrix of the pseudo spin at site 
i. The system starts the time evolution in a product state 


m 

l'I'o)= y] (23) 


with = cos{0j /2) and = sin(^j/2). We can formally write the 
time-dependent density matrix as 


m 

M j ' 

(24) 

Here, E[{a}] = hiai + . Vijaiaj. We can now trace out all 

but one site to obtain the reduced density matrix of a single pseudo-spin 
at site i in the basis of the eigenvalues of af : 


Pn 

= cos^ Oi/2 



(25) 

hi 

= sin^^i/2 

1 

Nf(t) 

(26) 

hi 

= ht = sin Oi 

1 

r7,= 1 

, (27) 

Nfit) 


where Un = E[\', {cr}] — E[],, {cr'}]. {a} and {a'} represent states of 
all other sites except site i. The number of frequencies involved (t) is 
roughly 2^^^\ i.e. all possible interactions between the l{t) spins that are 
significantly entangled with the observed spin at time t. More precisely, 
the number of frequencies Nf{t) = ^ Mut/h) j^e^sures the 

size of the entangled Hilbert space at the observation time. 

To find the time dependence of the CDW we have to transform back 
from the basis of af eigenvalues to the basis of af. This is achieved 
with a rotation around the axis by Oi 


(cr^(t)) = Tr(p(t)cr^) = cos^ Oi + sin^ Oi cos{uJnt). 

Nf{t) ^ 

(28) 

Hence, the fluctuations of the local imbalance between an even site and 
neighboring odd site behaves as 


ScTrmsit) 


e 


-hs(t) 



(29) 


while the fluctuation of the global imbalance (CDW order) is further 
suppressed by a factor of 1/ vT 


SXrmsit) 


—(-V 

^/zyutJ 


(30) 


We see that the decay of the fluctuations is intimately connected to the 
Int/r growth of the entanglement entropy, as it also ‘measures’ the 
number of spins coupled through the interactions after time t. 
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Figure 13: Decay of imbalance oscillations and entanglement entropy 
growth. A shows the time evolution of the imbalance with decreasing oscillation 
amplitudes for different interaction strengths for a system with 30% doublons, 
exhibiting a power law decay. B shows the connection of the decay exponents of 
the imbalance oscillation amplitudes and the slope of the logarithmic growth of 
entanglement entropy as a function of interaction strength. 

The above analysis is of a simplified effective model. However, 
the main conclusions are supported by direct simulation of the Hub¬ 
bard model on the quasi-periodic lattice using time dependent DMRG. 
Fig. pj] \ shows the temporal noise of the imbalance as a function of the 
time for different values of the interaction strength U/J. The fluctua¬ 
tions are measured by averaging them over a time window of T = 7 r 
around the time t. The results fit well to a power law decay. Fig. 
compares the fitted exponent s* to the slope of the logarithmic growth 
of the entanglement entropy showing the direct correspondence between 
the two. 

Thus, we conclude that measurements of the temporal fluctuations 
of the CDW order provide a viable experimental route to determine the 
bare localization length and distinguish the many-body localized state 
from an Anderson localized state of non-interacting particles. In the 
present experiment, this is, however, not possible, as the unavoidable 
averaging over many tubes suppresses the fluctuations below the detec¬ 
tion limit after only few oscillations. For future experiments, a single 
tube, or even single-site resolution would be desirable to overcome this 
limitation. We also note that the temporal fluctuations of the expectation 
value are different from shot-to-shot fluctuations at a given time which 
reflect the quantum uncertainty of the observable and would be finite 
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even in the infinite-time limit. 
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